1 



Progress of Theoretical Physics, Vol. 120, No. 5, November 2008 



Formation and Radiative Feedback of First Objects and First 

Galaxies 

Masayuki Umemura^'*\ Hajime Susa^', Kenji Hasegawa^', Tamon Suwa^', and 

Benoit Semelin^'^ 

^ Center for Computational Sciences, University of Tsukuha, Tsukuha 305-8577, 

Japan 

^ Department of Physics, Konan University, Okamoto, Kobe, Hyogo 658-8501, 

Japan 

^ Fujitsu Limited, Kawasaki, Kanagawa 211-8588, Japan 
^ LERMA, Observatoire de Paris, CNRS, 61 Av. de V Observatoire, 75014 Paris, 

France 

^ Universite Pierre et Marie Curie, 4 place Jussieu, 75005 Paris, France 

(Received ; revised ) 

First, the formation of first objects driven by dark matter is revisited by high-resolution 
hydrodynamic simulations. It is revealed that dark matter haloes of ^ lO^M© can produce 
first luminous objects with the aid of dark matter cusps. Therefore, the mass of first objects 
is smaller by roughly two orders of magnitude than in the previous prediction. This implies 
' that the number of Pop III stars formed in the early universe could be significantly larger 

than hitherto thought. Secondly, the feedback by photo- ionization and photo-dissociation 
photons in the first objects is explored with radiation hydrodynamic simulations, and it is 
demonstrated that multiple stars can form in a 10^ Mq halo. Thirdly, the fragmentation of 
an accretion disk around a primordial protostar is explored with photo-dissociation feedback. 
As a result, it is found that the photo-dissociation can reduce the mass accretion rate onto 
protostars. Also, protostars as small as O.SMq may be ejected and evolve with keeping their 
mass, which might be detected as "real first stars" in the Galactic halo. Finally, state-of-the- 
art radiation hydrodynamic simulations are performed to investigate the internal ionization 
of first galaxies and the escape of ionizing photons. We find that UV feedback by forming 
massive stars enhances the escape fraction even in a halo as massive as > 6 x 10^ M©, while 
it reduces the star formation rate significantly. This may have a momentous impact on the 
cosmic reionization. 

§1. Introduction 

The first generation objects and galaxies in the early universe are of great signif- 
icance as the progenitors of present-day galaxies, the sources for cosmic reionization, 
the origin of heavy elements, and the generators of seed black holes for supermassive 
black holes inhabiting galactic centers. 

The formation of first generation objects has been explored by many authors. 
Originally, the minimum halo mass that can undergo gravitational instability was 
estimated to be from 5 x IO^Mq at redshift z ^ 100 to 5 X lO'^M© at z ^ 10 
from analytic argument s,'^ where Mq is the solar mass. Later on, three-dimensional 
hydrodynamic simulation^ showed the minimum mass to be from lO^M© at redshift 
Z 100 to 10*^ M© at z ^ 10. Then, full cosmological hydrodynamic simulations have 
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been performed to show that the minimum mass is a weak function of redshift, which 

is ^ 7 X 10^ Mq in halo mass and ^ lO^M© in baryonic massP As for the secondary 

cohapse of first objects under the influence of preforming stars, the radiative feedback 

through photo-ionization of hydrogen, photo-dissociation of hydrogen molecules, and 

their combination are important. This issue has been explored by the radiation 

hydrodynamics, where hydrodynamic simulations are coupled with radiative transfer 
calculationsP'EJ'IS 

Also, in the last decade, the mass of Population III (Pop III) stars has been 
studied extensively. It has been shown that first stars are likely to be as massive as 
100 - 1OOOM0.[IS1)'IID'[I3'EI) On the other hand, if the collapse of a first object into a 
flat disk is considered, first stars may form in a bimodal fashion with peaks of ^ IMq 
and several lOOM©.^^^ Besides, the mass of massive Pop III stars could be reduced to 
20 — 4OM0 by cosmic variance, external feedback, or internal feedback. These 
works focused on the runaway collapse phase of first objects. Very recently it has 
been shown that in the accretion phase the disk around a first protostar can fragment 
into smaller pieces, eventually allowing the formation of less massive stars down to 
^ IMq or subsolar value .'I^'flSl) 'ESI) Even in such a fragmented primordial disk, stars 
that form earlier can make a significant impact through radiation hydrodynamic 
feedback on subsequent star formation. 

The first generation and second generation stars are assembled into first galaxies. 
Ultraviolet (UV) radiation from massive stars in first galaxies can play an important 
role not only on the following star formation history but also as the source for 
cosmic reionization. The escape of ionizing photons from first galaxies is regulated 
by the ionization of interstellar medium. This issue should be also investigated with 
radiation hydrodynamics. 

In this paper, we present recent high-resolution hydrodynamic simulations on 
first objects and radiation hydrodynamic simulations on the radiative feedback in 
these objects. Also, the fragmentation of an accretion disk around a first protostar is 
explored with photo-dissociation feedback, and radiation hydrodynamic simulations 
are performed to investigate the internal ionization of first galaxies and the escape of 
ionizing photons. In §2, a dedicated simulator for radiation hydrodynamics, which 
has been developed for the present purpose, is described. In §3, the collapse of 
first objects driven by dark matter cusps is explored. The radiation hydrodynamic 
feedback there is studied in In §4. In §5, the fragmentation of a first protostellar 
disk and radiative feedback by photo-dissociation photons is explored with radia- 
tion hydrodynamic simulations. In §6, the photo-ionization of first galaxies and the 
escape fraction of ionizing photons are analyzed based on radiation hydrodynamic 
simulations. §7 is devoted to the summary. 

§2. FIRST Project 

In order to explore extensively the hydrodynamic and radiation hydrodynamic 
processes during the formation of first objects and first galaxies, we have built up a 
dedicated simulator called FIRST, on the basis of the F/i?5T project.^^^ 

The FIRST project was initiated by a Specially Promoted Research in Grants- 
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Fig. 1. The F/i?5'T Simulator composed of 256 nodes [left) and an accelerator, Blade-GRAPE, for 
gravity calculations, which is embedded in each node {right). 

in- Aid for Scientific Research by MEXT over four years (2004-2007) with a budget 
of JPY 329.5 mihion (US $4.1 mihion), and has been continued with a Grant-in- Aid 
for Scientific Research (S) (2008-2012) by JSPS with a budget of JPY 73.1 mihion 
(US $0.9 mihion). In this project, we constructed a new hybrid simulator FIRST 
(Fusional Integrator for Radiation- hydrodynamic Systems in Tsukuba University). 
The F/iJS'T simulator is a large-scale hybrid PC cluster, where each node possesses a 
newly-developed board for gravity calculations, Blade-GRAPE. The Blade-GRAPE 
is composed of four GRAPE-6A chips.^^ The theoretical peak performance of one 
Blade-GRAPE board is 136.8 Gflops. Each board has 16MB of memory and can 
treat 260,000 particles simultaneously. The Blade-GRAPE is directly connected via 
PCI-X bus. Each server PC is equipped with a multi-port Gigabit Ethernet NIC 
that is connected to a special interconnection network using commodity Ethernet 
switches. 

The first version of Blade-GRAPE works with 32 bit and 33 MHz (PCI). Then, 
it has been improved to a 64 bit and 100 MHz version (PCI-X), Blade-GRAPE X64. 
Using Blade-GRAPEs, we have constructed a hybrid PC cluster system composed 
of 256 nodes, the FIRST simulator (Fig. [2j). The system possesses 224 Blade- 
GRAPE X64 boards and 16 Blade-GRAPE boards. The host PC cluster node is 
a 2U-size server PC (HP ProLiant DL380 G4) that has dual Xeon processors in a 
SMP configuration. As a result, the total performance of the FIRST simulator is 
36.1 Tflops, where the host PC cluster is 3.1 Tflops and the Blade-GRAPEs are 
33 Tflops. All nodes are connected uniformly to each other via a multi-port Gbit 
ethernet interconnect switch. The total memory of the FIRST simulator is 1.6 
TB. Also, the Gfarm Grid file system ( |http :/ /datafarm. apgrid.org/index. en. html) 
is installed. With Gfarm, a total storage of 89.2 TB is available as a seamless file 
server. 

In this paper, we present hydrodynamic and radiation hydrodynamic simulations 
on the formation of first objects and first galaxies with the FIRST simulator. 
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§3. Collapse of First Objects by Dark Matter Cusps 

The efficiency of the formation of first objects is extremely significant in that they 
are responsible for the reionization history in the universe and the cosmic chemical 
enrichment. Yoshida et al. (2003)^^ have shown that the minimum mass is a weak 
function of redshifts, which is ^ 7 X IO^Mq in halo mass and ^ lO^M© in baryonic 
mass. The minimum mass is basically determined by the thermal instability via 
hydrogen molecules (H2). H2 molecules form in the non-equilibrium processes with 
the catalysis of free electrons as 

e~ + H R- + hu 

H- + H ^ H2 + e-. 

If the temperature exceeds lO^K around the number density of 1 cm~^, the formation 
rate of H2 molecules is increased, resulting in the H2 abundance of ^ 10~^. 
Then, thermal instability occurs and the temperature decreases down to ^ 300K 
with increasing density. Then, the runaway collapse of the cloud proceeds through 
gravitational instability. The simulations by Yoshida et al. (2003) demonstrated 
that for the halo mass below 7 x lO^M©, the virial temperature cannot be raised 
up to ^ lO^K, and therefore the thermal instability by H2 cooling does not occur. 
The mass resolution of their simulations is ^ 3OM0. Recently, Umemura et al.'^''^ 
have found by simulations with higer mass resolution that significantly smaller mass 
haloes allow the collapse of primodial clouds via the H2 cooling instability. The 
major difference from previous simulations is the fact that the mass resolution is 
much higher and not changed throughout the simulations, in contrast to the adaptive 
change of resolutions in previous works. Umemura et al. argued that the runaway in 
a smaller dark halo could be induced by the cusp of dark halo. In other words, the 
resolved cusp potential can raise the central temperature of cloud up to - lO^K and 
leads to the thermal instability via H2. In this paper, to verify this conjecture, we 
compare the numerical results with different mass resolution, and analyze the dark 
matter cusp-induced collapse. 

3.1. Set-up of Simulations 

The formation of first objects from primordial gas is driven by dark matter fluc- 
tuations and cooling by H2 molecules. In order to investigate the dependence on 
dark matter potential, we perform high-resolution cosmological hydrodynamic sim- 
ulations in a standard ^CDM cosmology with the parameters of {Oa, ^m, ^h^h) = 



Table L 


Models and 


mass resolution 


Model 


(ArDM,A^b) 


(mDM,mb) 
[Mo] 


R64 


(64^64^) 


(24.1,5.08) 


R128 


(128^128^) 


(3.01,0.64) 


R256 


(256^2562) 
(512^5122) 


(0.38,0.079) 


R512 


(0.047, 0.0099) 



(3-1) 
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Fig. 2. Cumulative mass of dark matter and baryons is shown against radius at z — 16 for different 
mass resolution simulations. The left panel shows model R64, the middle panel does model 
R128, and the right panel does model R256. A straight line shows the mass corresponding to 
the density distribution of p oc r~^'^. The Jeans mass including dark matter that is defined by 
()3 4p is also shown. 

(0.72,0.24,0.042,0.73). First, we perform pure A^-body simulations for dark matter 
in a larger box with a commoving linear scale of IMpc and find the highest den- 
sity domain of (60kpc)^. For this domain, the evolution of density fluctuations is 
recalculated including baryons from the recombination epoch {z = 10^). We use a 
Particle-Particle-Particle-Mesh (P^M) scheme for gravity force calculations, and the 
baryonic component is treated with the Smoothed Particle Hydrodynamics (SPH) 
method. Radiative heating/cooling and chemical reactions for e~, H, H+, H~, H2, 
H2" are included, where helium and deuterium reactions are dismissed. The direct 
part of self-gravity is calculated with Blade-GRAPE in the FIRST simulator. The 
simulations do not use any zoom-in technique, but the mass resolution is uniform 
over the whole computational domain and throughout the evolution. The system 
mass is 6.3 x lO^M© in dark matter and 1.6 x 10^ Mq in baryonic matter. The 
models with different resolutions are summarized in Table [B 

3.2. Growth of dark matter cusps 

A dark matter fluctuation in the computational domain enters a non-linear stage 
around z = 30. In the course of the non-linear evolution, dark matter forms a 
virialized halo, and simultaneously a central cusp develops. The growth of the cusp 
depends on the mass resolution. In Fig. 13. 2^ the mass distributions are shown for 
different mass resolution simulations. In model R256, the cusp intrudes to radii of 
^ Ipc. The density distribution is well fitted by p oc r~^'^. In a low or intermediate 
resolution case (model R64 or R128), the cusp is smoothed out at innermost regions. 
The dark matter density profile with a cusp is often expressed by 

p{r) = Psir/r.r'^il + {r/r,)'f--y\ (3-2) 

wehre is a characteristic inner radius and ps is the corresponding inner density. 
In the NFW profile, (tt, 6) = (1? 1) and the mass enclosed within r is M(r) = 
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Fig. 3. Evolution of baryonic density fluctuations from z — 90 to z — 17 for model R256. Colors 
show the density according to the legend of colored bar. 

47rr^ln[l + (r/r^)] - {r/rs)/[l + (r/r^)] or in the Moore's profil^ (a, 6) = (1.5, 1.5) 
and the mass is 

M(r) = (87r/3)p,rfln[l + (r/r,)i-^]. (3-3) 

The pure A^-body simulations of dark matter hitherto performed have shown that 
the cusp profiles are between the NFW profile and the Moore's one.E^'l^'EBSa 
The profiles of cusps {p oc r~^'^) in the present simulations are close to the Moore's 
profile, although the preset results are slightly steeper. (This may be the effect of 
dissipation of baryonic component.) For a virialized halo with 10^ M©, the cusp 
mass is ^ 10^ Mq if the concentration parameter is c = Vyij-frg ^10 with the virial 
radius ryir- If the mass of dark matter particle is higher than Mq, a cusp with 
^ W^Mq is made of less than several hundred particles. The two-body relaxation 
proceeds in a timescale of O.lA^tdyn/lnA^, where N is the number of particles and t^yn 
is the dynamical time. If < 10^, the two-body relaxation timescale is less than 
Stdyn- Thus, the dark matter cusp is inevitably erased by the two-body relaxation. 
Actually, in model R64 {ttidm = 24.IM0) or R128 {ttidm = S.OIM©), the mass 
distribution decreases steeply toward the center. This meams that the central cusps 
are smoothed out by the two-body relaxation. On the other hand, in model R256 
{t^dm — O.38M0), a cusp grows down to ^Ipc. Consequently, the central dark 
matter potential becomes deeper. 
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Fig. 4. Spatial density distributions of baryons and dark matter z — 12 are shown for different 
resolution simulations. The left panel shows model R64, the middle panel does model R128, and 
the right panel does model R256. Colors shows baryonic density following the attached legend, 
and contours show the levels of dark matter density. 



In Fig. \i.2\ we also show the Jeans mass defined by 

where Cg is the local sound velocity, and Mdm and Mb are the cumulative mass of 
dark matter and baryons, respectively. As seen clearly, for higher mass resolution 
the Jeans mass is reduced in innermost regions owing to the dark matter cusp. The 
central dark matter potential is responsible for the temperature, which is directly 
related to the thermal instability by H2 molecules. 

3.3. Thermal instability induced by a dark matter cusp 

In Fig. El the evolution of baryonic density fluctuations induced by dark matter 
is shown from z = 90 to z = 17 for model R256. If the dark matter particle mass 
is less than IM0, the cusp with the density profile of r~^'^ is resolved down to 
< 10^ Mq. In contrast to the previous works, multiple peaks can collapse induced by 
dark matter cusps. In Fig. [U the density distributions at z = 12 are compared among 
different resolutions. If the dark matter resolution is lower than ttidm = lOM©, the 
cusp cannot be resolved sufficiently in a halo of 10^ M©, and small-scale peaks readily 
merge into a larger peak, and consequently only one baryon condensation collapses 
within a halo. This evolution is basically equivalent to the previous works. But, in 
model R256 (moM = O.38M0), a double peak forms, where the mass of each peak 
is ^ IO^Mq in dark matter and ^ IO^Mq in baryons. The separation of peaks is 
roughly 60pc. Unless they undergo the thermal instability, the collapse of each peak 
bounces and the peaks may eventually merge into one larger peak. Using eq fj3-3IK 
the baryonic gas temperature in the cusp region can be evaluated as 

_ ^ ^ ( P. \ (-r^y(-^y\ ,3.5) 

^ r kB \1.6xl0 ^^gcm-^y \10pc/ \10pc/ 
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where /i is the mean molecular weight, rrip is the proton mass, and A:^ is the Boltz- 
mann constant. Thus, if the cusp potential is resolved, the temperature can be 
raised up to > lO^K, and then thermal instability through H2 cooling can occur. In 
Fig. [5l the temperature and H2 fractions in model R256 are shown versus baryon 
number density. Top panels present those for all particles, middle panels for the 
particles in the highest peak, and bottom panels for the particles in the 2nd highest 
peak. This figure shows that H2 fractions reach a level of several 10""^ around the 
density of nn ^ lcm~^ in the first and 2nd peaks. Then, the temperature decreases 
via thermal instability in both peaks. What is important is that although the halo 
mass of peaks is of the order of ^ lO^M©, the temperature reaches over ~ lO^K 
around nn ^ lcm~^ owing to the deep potential of dark matter cusps. Eventually, 
dark haloes of ^ 10"^ Mq can produce primordial objects. It means that the mass of 
first generation objects can be smaller by roughly two order of magnitude than in 
previous prediction. Hence, we conclude that the number of Pop III stars formed in 
an early universe can be significantly larger than in the previous prediction. 

We have demonstrated that the resolution of dark matter cusps is crucial for 
the thermal instability. Here, we make the convergence test regarding the resolution 
to verify the present results. We perform a ultra-high resolution simulation (model 
R512). In Fig. El we compare the results between models R256 and R512. We 
confirm that both give basically the same results. This means that the dark matter 
cusp resolved down to Ipc is important for thermal instability. The further inner 
cusp is not responsible for a significant increase of temperature. 

§4. Radiation Hydrodynamic Feedback in First Objects 

In the previous section, we have shown that two peaks with separation of 60pc 
collapse independently. Although each peak could form a fairly massive star of 
> 4OM0,i^) there is a difference of free-fall time of ^ 5 x lO^yr. Hence, the first peak 
forms a massive Pop HI appreciably ealier than the 2nd peak, and then the first star 
irradiates the second collapsing peak with strong UV radiation. Thus, whether the 
2nd peak successfully forms a star depends upon the strength of radiative feedback 
from the first star. The physical mechanism of radiation hydrodynamic feedback has 
been studied by Susa & Umemura (2006), Susa et al. (2009) ,13 and Hasegawa et 
al. (2009) .'^ We show below the basic physics of radiation hydrodynamic feedback 
and explore the survival of the 2nd peak under the feedback. 

4.1. Basic Physics 

The properties of Pop HI stars in the range of 2^Mq < < I2OM0 are 
summarized in Table HH which are taken from Schaerer (2002).^^^ Since Pop HI stars 
emit strong ionizing and photo-dissociating radiation, the radiation hydrodynamic 
feedback is regulated by the propagation of an ionizing front and the shielding from 
photo-dissociating radiation. If a collapsing core is irradiated by an ionizing source 
located at a distance the propagation speed of the ionization- front (I- front) in 
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Fig. 5. Temperature (left) and H2 fractions (right) versus barynon number density in model R256 
are shown for all particles {top panels), the particles in the highest peak {middle panels), and 
the particles in the second highest peak {bottom panels). 



the core is given by 




if recombination is neglected, where A^ion is the number of ionizing photons per unit 
time and ricore is the number density in the cloud core. The sound speed in regions 
cooled by H2 IS (1\ ^ 1km s ^, while that in the ionized regions is fl2 ^ 10km s-\ 
If the density of the cloud core is low or the ionizing radiation flux is strong, then 
> 2a2 and therefore the I-front becomes R-type. If we focus on a core collapsing 
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Fig. 6. Baryonic density distributions are compared at z = 15.6 between model R256 (left) and 
R512 (right). Colors shows baryonic density following the attached legend. 

in a self-similar fashion, the core size r^ore is on the order of ttitff, where tff (— 
core) is the free-fall time. Then, the propagation time of R-type front across 
the core satisfies tip = r^ore 

/viF < (ai/2a2)tff < tff. This means that an R-type 
front sweeps the core before the core collapses in the free-fall time. Thus, the core is 
likely to undergo photo-evaporation. On the other hand, if the density of the cloud 
core is high enough or the source distance is large, then v^^ < a\/2a2 and a D-type 
I-front emerges. The propagation time of a D-type front across the core satisfies 
tiY > (2a2/ai)tff > tff. Thus, the core can collapse before the I-front sweeps the core. 

However, the above arguments are based on the assumption that the ionizing 
photon flux does not change during the propagation of I-front. The core could be self- 
shielded from the ionizing radiation if 7V(7rr^Qj.g/47ri?^) < 47rr^Qj.gn^Qj.gaB/3, where 
Q^B is the recombination coefficient to all excited levels of hydrogen. The critical 
density for self-shielding is given by 

/ 37Vion [G^Y^ 

^shield - T^—^ V 



Table II. Properties of Pop III stars 



Mass 


Teff [K] 




ivLw[erg/s] 


I2OM0 


9.57 X 10^ 


1.069 X 10^° 


5.34 X 10^^ 


8OM0 


9.33 X 10^ 


5.938 X 10^^ 


3.05 X 10^^ 


4OM0 


7.94 X 10^ 


1.873 X lO'^^ 


1.17 X 10^^ 


25M0 


7.08 X 10^ 


5.446 X 10^^ 


3.94 X 10^^ 



a) Effective temperature 

b) Number of ionizing photons emitted per second 

c) Luminosity of photo-dissociation radiation at Lyman- Werner (LW) band 
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Fig. 7. Time variations of gas density, temperature, H2 fractions, and ionization degree along the 
line connecting two peaks are shown at the epochs from O.lSMyr to 2.5Myr. The distance is 
measured from the first peak. The 2nd peak is located at ^ 60pc 



If ricore > ^shield 5 the ionizing photon flux diminishes signiflcantly during the I- front 
propagation. Hence, even if the I-front is R-type on the surface of the cloud core, 
the front can change to M-type accompanied by a shock, and eventuahy to D-type 
inside the core.l^ 

In contrast to ionizing radiation, H2 dissociating radiation in LW band (11.26- 
13.6 eV) is less shielded. The self-shielding of LW band flux (Flw) is expressed 
by 

i^LW = i^LW,0/sh(A^H2,14) (4-3) 

where Flw,o is the incident flux, A^H2,14 
column density, and 

Hence, if the column density of H2 molecules ahead of I-front is high enough, H2 
dissociating radiation can be shielded. Since H2 molecules form with the catalysis of 
free electrons, the shielding of dissociating radiation is intimately coupled with the 
propagation of I-front. 



A^H2/10 cm is the normalized H2 
" - \ (4-4) 

X > 1 ^ ^ 
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4.2. START: Accelerated ratiation hydrodynamic scheme 

In order to perform the radiation hydrodynamic (RHD) simulations more ef- 
fectively, recently we have developed a novel RHD code called START (SPH with 
Tree based Accelerated Radiation Transfer), which is designed to solve the transfer 
of UV photons from numerous sources. In START^ the optical depth between a 
radiation source and an SPH particle is evaluated with the same method as in the 
RSPH code's^ In contrast to the corresponding method proposed by Kessel-Deynet 
& Burkert (2000) ,^3 the evaluation of the optical depth between a radiation source 
and a target SPH particle is performed only once by using information of its up- 
stream particle. Thus, the cost for each ray-tracing turns out to be proportional to 
^SPH, where A/'sph is the numbers of SPH particles. In addition to the reduction of 
the cost for each ray-tracing, the effective number of radiation sources for each tar- 
get SPH particle is diminished by utilizing an oct-tree structure of radiation sources. 
Similar to the Barnes-Hut tree method^^^ that is frequently used for calculating the 
gravitational force, if a cell containing some radiation sources is far enough from a 
target SPH particle, all radiation sources in the cell is regarded as a virtual lumi- 
nous source. As a result, the computational cost to evaluate optical depths from all 
radiation sources to all SPH particles is proportional to 0(A^sph log A^s). With this 
code, we can precisely solve the transfer of recombination photons, which is often 
treated with the on-the-spot approximation. Another important and useful feature 
of START is that the spatial resolution is adaptively enhanced in high density re- 
gions, since SPH particles are directly used for the radiative transfer grids. Here, 
using START, we consistently solve the gas and dark matter dynamics coupled with 
the radiative transfer of UV photons and non-equilibrium chemical reactions for 6 
species: e~, H+, H, H~, H2, and H2". 

4.3. Radiation hydrodynamic feedback 

We use the results obtained in the previous section, where two dark matter 
haloes host gas clumps, and the separation between the halos is ^60pc. The whole 
simulation box size is 2kpc (physical) on a side. The particle mass is O.OSMq for an 
SPH particle, and 0.38Mq for a DM particle, respectively. We assume a source Pop 
HI star with a mass of I2OM0, and place it at the highest density peak. Then, we 
start the RHD simulations, where not only direct UV photons from the source Pop III 
star but also diffuse photons produced via recombination processes are considered. 
After the lifetime of the source star (t = 2.5Myr), direct UV photons from the source 
star is turned off, but the transfer of diffuse photons is continuously solved. After 
the source star dies, the supernova explosion may make an impact on the collapsing 
core, which is not treated in this paper and will be explored elsewhere. 

The evolution from O.lSMyr to 2.5Myr obtained by the radiation hydrodynamic 
simulation is presented in Fig. [71 where time variations of gas density, temperature, 
H2 fractions, and ionization degree along the line connecting two peaks are shown. 
The distance is measured from the first peak, and the 2nd peak is originally located 
at ^ 60pc. In the early stage before t = 0.75Myr, a shock precedes a D-type I-front. 
Before the irradiation of UV, the H2 fraction at the 2nd peak is several 10~^, as 
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Fig. 8. Snapshots for the gas density {upper row), and the ionization degree {lower row) at 2.5Myr. 
The initial distribution at t = O.OMyr corresponds to the rightmost panel of Fig. [4l The origin, 
(x,y) = (0,0) [pc], is set to be the position where the highest peak was originally located, and 
the 2nd peak was originally located at (x,y) = (20,50) [pc]. The contours show the dark matter 
distributions. The left panels are the results under the on-the-spot approximation, while the 
right panels show the results of full radiation-hydrodynamic simulations solving the transfer of 
diffuse photons. 



shown in Fig. O However, as shown in Fig. [71 the the H2 fraction at the 2nd peak 
is reduced to a level of 10~^ due to photo-dissociating radiation, although the H2 
fraction near the I-front is slightly enhanced. After this stage, the I-front changes 
into R-type. The evolution until this early stage is very similar to the corresponding 
results shown by Kitayama et al. (2004)^^ or Yoshida et al. (2007).^^ However, 
unlike these previous studies, the I-front changes into D-type again at t ^ 1.25Myr 
due to the presence of the 2nd peak. Then, an H2 shell with an H2 fraction of ^ 10~^ 
forms ahead of the 2nd peak thanks to the catalysis of free electrons. This H2 shell 
shields photo-dissociating radiation significantly, and consequently the H2 fraction 
at the 2nd peak is restored to ^ 10~^. Simultaneously, the envelope of the 2nd peak 
is stripped by the shock associated with the I-front. 

The ionization around the 2nd peak is dependent on the treatment of UV ra- 
diation transfer. In Fig. [SI the results under the on-the-spot approximation are 
compared to those of full RHD simulations solving the transfer of diffuse photons at 
the epoch of 2.5Myr. As shown in the top panels of Fig. [HI the diffuse UV radiation 
hardly affects the gas distribution in the 2nd peak, because the mean free path of 
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Fig. 9. The late-phase evolution of the 2nd peak is shown at 5.5Myr, 32.5Myr, and 57.5Myr. The 
origin, (x,y) = (0,0) [pc], is set to be the position where the highest peak was originally located, 
and the 2nd peak was originally located at (x,y) = (20,50) [pc]. The gas density is shown in the 
upper row, and the ionization degree is shown in the lower row. 



ionizing photons is quite short there. Thus, the impact of the diffuse UV radiation 
on the 2nd peak is dimunitive. On the other hand, the difference of the ionized 
fractions in the outer envelope is obvious behind the 2nd peak. The fuh RHD simu- 
lations show that diffuse UV radiation can ionize the gas behind the 2nd peak, since 
the photon mean free path is relatively long in this region. Therefore, solving the 
transfer of diffuse UV radiation is of great significance in computing the ionization 
of the low-density gas, e.g., intergalactic medium. 

The later evolution of the 2nd peak is determined by the interplay of the shielding 
from photo-dissociating radiation and the stripping by an ionization-front shock. In 
Fig. [9l the later evolution in full RHD simulations is presented at 5.5Myr, 32.5Myr, 
and 57.5Myr, where the distributions of density and ionization degree are shown. 
After t =2.5Myr, no UV radiation is emitted from the source. Thus, the recom- 
bination gradually increases the neutral fractions. But, diffuse UV radiation from 
the recombination ionizes the shadowed region behind the 2nd peak. Although the 
central density of the 2nd peak decreases via the phot oionizat ion and a shock asso- 
ciated with the ionization-front, the 2nd peak collapses eventually. In Fig. (TUl the 
the density of the 2nd peak is shown against temperature. This figure shows that 
the 2nd peak continues to collapse after the thermal instability, irrespective of the 
radiative feedback. 
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Fig. 10. Temprature versus density at the end of the simulation. 



4.4. Analytic Criteria 

Here, we derive analytic criteria for the feedback, based on the argument by 
Hasegawa et al. (2009).^^ Susa (2007) ^"^^ has explored the photo-dissociation feed- 
back of a Pop III star with I2OM0 on a neighboring prestellar core by RHD simu- 
lations without ionizing radiation. Susa (2007) has found that a condition for the 
collapse of a neighboring core is approximately determined by t^is = ^fr, where t^is is 
the photo-dissociation timescale in the core and tff is the free-fall timescale. Using 
this condition, the critical distance i?cr,d5 above which a neighboring core is able to 
collapse, is given by 



where Llw? ^c? and are the LW luminosity of the source star, the number density 
of the core, and the temperature of the core, respectively. This equation involves the 

self-shielding effect by the core. The dependence on the core temperature basically 

1/2 

originates in the core radius (oc Tc ) and a H2 formation rate in the core (oc Tc). 
Hence, the self-shielding effect is weaker for a lower core temperature. The boundary 
between the collapse irrespective of ionizing radiation and that with the aid of ion- 
izing radiation can be determined by i?cr,d- Iii addition, i?cr,d gives a good estimate 
for less massive source star cases, because the ionizing radiation is relatively weak 
for less massive stars. 

If ionizing radiation is included, we should also incorporate the shielding effect 
by an H2 shell, which is formed with the catalysis of free electrons. Here, we derive 
a new criterion including this effect. Since a cloud collapses in a self-similar fashion 
before UV irradiation, the density profile of the outer envelope in the cloud is well 
expressed as 




(4-5) 




(4-6) 
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where is the core radius which roughly corresponds to the Jeans scale; 



where denotes the proton mass. Assuming that the thickness of the H2 shell is 
determined by the amount of ionized gas in the envelope and the H2 fraction in the 
shell is constant, the H2 column density of the shell A^H2,sh is given by 

Jd J^J^sh 

where i?sh and yH2,sh are the distance between the cloud core and the H2 shell, and 
the H2 fraction in the shell, respectively. Here, D^y^ is set to be the position where the 
number of recombination per unit time in the ionized region around a source star 
balances the number rate of incident ionizing photons, since the H2 shell appears 
ahead of the ionization front. Hence, D^h satisfies 



— 27rQ;B / n{r)'^r^dr 
Jd 

) - A 



27r«BnM ^r.rf''' - (4-9) 



Using equation fl4-8l) and fl4-9l) . we obtain 



12 2 f N- \^ 
Nn„s^^yn„s^niriD-^[-^\ . (4-10) 

Because of the intense LW radiation, the H2 abundance at the shell is in chemical 
equilibrium. Therefore, yH2,sh is given by 

yH2,sh = 7 , (4-11) 

^dis 

where ye,sh is the electron fraction at the H2 shell and /ch- is the reaction rate of the 
H~ process by fj31l) . In this case, we should consider the self-shielding effect by the 
shell itself. As a result, these rates are 

= 1.0 X 10"^^rshcm"^s"\ (4-12) 

fcdis = 1.13 X 10^FLw,sh/s [^i^) (4-13) 

where Tgh and FLW,sh are the temperature at the shell, and the LW flux from the 
star in the absence of shielding effect, FLW,sh — /Ati{D — i^sh)^- /s is the self- 
shielding function, which is given by (j4-4p . Combining equations (I4-I0p - (I4-I3p with 
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assumption of ye,sh — 0.05 and Tgh = 2000K as shown in the present numerical 
results, we have 



2/H2,sh = 1.0 X lo- 



tion \ / ^LW 



lO^Os-i / \5x 1023ergs-i 



1 2 



Tc y3 f D Y 

300Ky V40pcy ' 



(4-14) 




(4-15) 

It should be noted that A^H2,sh is independent of the core temperature Tc, but is 
determined solely by the ratio of TVion to Llw with strong dependence. 
Multiplying Llw in equation (|4-5I) by 



we obtain the critical distance as 
^cr,sh = 147pc 



5 X lO^^erg s~^y VlO^cm 



i 7_ 
16 



3 



300K 



(4-17) 



in which both shielding effects by the core and the H2 shell are taken into account. 
In particular, if A^H2,sh is larger than lO^'^cm"^, the critical distance can be expressed 
as 



^cr,sh = 78.8pc 



LW \ / ^ion 



5 X 1023erg s-i J I lO^^s 



103cm-3/ V300K 



(4-18) 



Based upon equation (14- 151) . the shielding effect by the shell becomes weaker as 
^ion/^LW decreases. Equation (|4-18l) gives a criterion of distance, above which a 
cloud core can collapse through the shielding of photo-dissociation radiation by an H2 
shell. The present results of radiation hydrodynamic simulations show that a cloud 
located at 60pc collapses under the UV feedback. This is broadly consistent with 
this analytic criterion. In the present simulations, the peak separation was ^60pc. 
Whether this separation is a typical value has not been revealed in the present 
analysis. It needs the case study using more realizations of random fluctuations, 
which is left in the future work. 
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§5. Fragmentation of a First Protostellar Disk and Radiative Feedback 

As stated in the Introduction, theoretical efforts in the last decade revealed 
that very first stars are more massive than stars forming in present-day galaxies. A 
few years before, it was believed that they are very massive (> lOOM©). However, 
recent numerical studies using the sink particle technique suggest that owing to the 
fragmentation in an accretion phase, they could be 0(10) xMq or less, and eventually 
they form a multiple stellar systeml^^'l^'l^ In some cases, it might be also possible 
for sub-solar mass first stars to form from such fragments. In addition, recent two 
dimensional radiation hydrodynamic simulations revealed that UV radiative feedback 
from first protostars is quite important to quench the mass accretion onto the first 
protostar.-*^^^ 

In this framework, we incorporate the UV radiative feedback effects into the 
three-dimensional chemo-hydrodynamical simulations, in order to follow the frag- 
mentation of an accreting gas disk in realistic circumstances. 



Radiative Feedback No Radiation 



Face-on View 




4000AU 



Edge-on View 



Fig. 11. Snapshots at 5500yrs after the first sink formation are shown. Red dots represent SPH 
particles, while black dots are sink particles. Left panels (face-on/edge-on view) are the results 
with radiative feedback, while the right panels without radiative feedback. 
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Fig. 12. Bird-eye view of the gas disk around the first sink at two epochs. Transparency describes 
the gas density, while the color shows the logarithmic fraction of H2 molecules (^/Hs)- At earlier 
epoch (1600yrs), conical zone at the center is highly dissociated (red gas), while the other region 
is still unaffected by the central source. On the other hand, at a later epoch (5500yrs), only 
dense cold streams (white/green gas) orbiting around the center sustain significant amount of 
H2, while the other less dense region is totally dissociated. 

5.1. Numerical method with sink particles 

We perform numerical experiments using the RSPH code developed by our- 
selves.133 We set the the initial condition of a cloud to be a Bonner-Evert sphere 
with the density of lO'^cm"^ and the temperature of 200K, at the "loitering" phase 
of a collapsing primordial prestellar cloud. We add rigid rotation with angular ve- 
locity of = 2 X 10~^^s~^, which is comparable to the rotation in cosmological 
simulations.^^ 

In order to follow the accretion phase, we introduce the sink particle technique 
in the RSPH code. In the case where the density at the position of an SPH particle 
exceeds lO^^cm"^, it is regarded as a sink particle. Also, particles that enter a sphere 
of radius 20AU centered on each sink particle are absorbed in the sink particle if they 
are gravitationally bounded. In the present simulations, the mass of the SPH particle 
is 6 X 1O~^M0. Following previous studies, the corresponding mass resolution is 
^ O.IM0, which is well below the Jeans mass of the disk at lOOOK, and that allows 
us to investigate the gravitational fragmentation of the disk. In the present work, 
we focus on the effects of Lyman- Werner radiation from the first protostar on the 
fragmentation of the disk, and the resultant accretion rate onto the first protostars. 
We solve non-equilibrium processes for H2 formation. The transfer of the Lyman- 
Werner radiation is solved utilizing the self-shielding function.'^ The luminosity and 
the effective temperature of the protostars are assessed by regarding the masses and 
the mass accretion rates of the sink particles as those of the protostars. In paticular, 
we interpolate the tables of luminosity and temperature generated by Hosokawa et 
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Fig. 13. Time evolution of total sink mass. Two curves correspond to the cases with/without 
radiative feedback. 

al.ISS Although the effects of photoionization are also important, we could not treat 
them due to the lack of spacial resolution. We leave them in future works. 

5.2. Radiative feedback in a first protostellar disk 

Fig. [TT] shows the numerical results, where the snapshots of distributions are 
presented at 5500yrs . Left panels show the edge-on/face-on views of the gas disk in 
the inner (4000 AU)^ box of the simulated region with radiative feedback, whereas 
the right panels for the case without radiative feedback. In both cases, we observe 
that several sink particles are formed (black circles), while the gas particle (red dots) 
distributions are more extended in a case with radiative feedback. This difference 
is a consequence of the H2 photo-dissociation process. Two panels in Fig. [12] show 
the bird-eye views of the gas disk around the protostars at two epochs (IGOOyrs 
and 5500yrs). At the later epoch (right panel), the polar region of the gas disk is 
highly dissociated by the radiation from the sink particles. The dissociated region 
is heated up to 7000 — 8000K, due to the accretion shocks and the chemical heating. 
In the absence of photo-dissociating radiation, H2 dissociation cooling should almost 
cancels the H2 formation heating. But, in the photo-dissociated region, H2 formation 
heating is of great significance, since H2 dissociation cooling is negligible in such 
environments. 

The "heated bubble" originating from the absence of coolant as well as the 
chemical heating by photo-dissociation reduces the mass accretion rate onto the sink 
particles. In fact, the total accreted mass after 24000 yrs is ^ 50% smaller than the 
mass without radiative feedback (Fig. [T3|). Thus, the photo-dissociation is one of the 
important effects to control the final mass of the first stars. On the other hand, it is 
worth noting that the photo-dissociation feedback cannot quench the fragmentation 
of the gas disks. We observe the fragmentation of the disk in both cases with and 
without radiative feedback, and the numbers of the sink particles are ^ 10. 

We also find that a few sinks are kicked away from the central part of the 
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host gas cloud via the gravitational three-body interaction. As a result, the mass 
accretion onto these sinks is basically shut off. We find that some of these sinks do 
not grow to more than O.SMq. The lifetime of main sequence of such low mass stars 
is ^ 2 X lO-'^^yr. If we regard these sink particles as "stars", we might be able to find 
the "real first stars" in the Galactic halo. 

§6. Evolution of First Galaxies by Internal UV feedback 

UV radiation affects the chemo-thermal and dynamical evolution of galaxies 
through the photo-ionization, -heating, and -dissociation processes. Although the 
importance of such UV feedback on structure formation in the universe have been 
pointed out already in 1980's, no calculation of three-dimensional hydrodynamics 
coupled with radiative transfer has ever been realized until 2000's. In the last 
decade, it became possible to perform three-dimensional radiation hydrodynamics 
(3D-RHD) simulations, thanks to the developments of computers and calculation 
algorithms.ES The START scheme enables us to solve 3D-RHD for numerous ra- 
diation sources.SD In this section, we present the results of RHD simulations on 
the ionization of first galaxies and the escape of ionizing photons, and thereby 
demonstrate the necessity of RHD simulations. Here, a ^CDM cosmology with 
{Ha, ^2m, ^b, h) = (0.73, 0.27, 0.049, 0.71) is employed. 

6.1. Determination of radiative reaction rates 

Once optical depths are determined, radiative reaction rates such as the photo- 
ionization and photo-heating rates can be evaluated. For instance, the ionization 
rate for z-th chemical species (i.e., HI, Hel, or Hell) at a position r is generally given 



where n^, 2/^,2, and ai^^i are the number density, the Lyman limit frequency, and the 
cross-section for i-th species. Ii^^q is the intensity of incident radiation at frequency 
z/, and Tjy is the optical depth. Here, we should carefully treat this equation, because 
it becomes numerically zero where riy{r) is much larger than unity. It implies that 
the ionization rates become zero even if the number of ionizing photons is enough 
to ionize the medium. This leads to the unphysical feature that the ionization front 
does not propagate into an optically thick media. To avoid this problem, the form 
given by 



is often used instead of equation (|6-ip so as to conserve the photon number. The 
photo-heating rates can be similarly evaluated by multiplying the term in the integral 
by {hiy - hvL^i). 



by 




(6-1) 




(6-2) 
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6.2. Chemo-Thermal Evolution coupled with RHD 

We update the information of radiation field at every period of tjiT = 0.1 x 
min(trec5 ^hydro) by solving radiative transfer (RT), where tree and thydro 
spectively the timescales for the recombination tree — l/na(T) and hydrodynamics 
^hydro — ^/( I ^1 + (T)) , whcrc n, T, /i, V, and Cs(T) are the number density, temper- 
ature, smoothing length, velocity, and sound velocity of a SPH particle, respectively. 
In a case of thydro > ^rec? we repeatedly solve the RT at an interval of tree during 
^hydro- Order to determine tree 5 we always use only particles on which the rela- 
tive change of ionized fraction during previous tRT is greater than 10 percent. This 
treatment allows us to effectively reduce the computational cost due to the following 
two reasons. First, it is often true that high density regions are well self-shielded, 
thus tree should be determined among SPH particles that reside outside of the self- 
shielded regions, otherwise we have to solve RT many times to no avail. Secondly, at 
a later phase of expansion of an HII region, which is called D-type ionization front, 
the ionization equilibrium is almost fulfilled in the HII region, and the change of 
radiation field is caused mainly by the hydrodynamic motion driven by the thermal 
pressure. 

When the UV feedback is included, the chemical timescale (tehem,i — '^i/'^i) 
and thermal one (tthermai — T/T) are generally shorter than tuT- In addition, 
the set of equations which determine the abundance of chemical species are stiff. 
Therefore, we implicitly solve the chemo-thermal evolution with subcycle timestep 
tsub — iiiiii(tehem 5 ^thermal)- Here notc that the ionization rates also vary during tRT^ 
since they are multiplied by ni — fin/ {jimYi)^ where fi and ji are the fraction of z-th 
species and the mean molecular weight, respectively. 

6.3. Feedback and escape of ionizing radiation in first galaxies 

It is important for revealing the cosmic reionization history to clarify how many 
ionizing photons can escape from galaxies. Previous numerical simulations that 
evaluate the escape fractions of ionizing photons from galaxies can be categorized 
into two types: One is RDH simulations in which UV feedback is consistently taken 
into account and the other is non-RHD simulations in which the transfer of 

ionizing photons is solved without the back-reactions on chemodynamics caused by 
UV feedback.SSJ'SB The UV feedback is expected to play a crucial role on the 
determination of the escape fraction, since the gas in a less massive galaxy is easily 
evacuated by the UV feedback. "^^^ Thus it is generally expected that the escape 
fraction from a less massive galaxy is enhanced by the UV feedback, although nobody 
has quantitatively clarified how much the UV feedback impacts the escape fraction. 
On the other hand, the escape fraction from a massive galaxy is often expected to be 
insensitive to the UV feedback, since the galaxy hardly loses its gaseous components. 
However, if the internal density structure of the galaxy is significantly changed by 
the UV feedback, the escape fraction is possibly changed. 

Here, we investigate the impact of UV feedback on the escape fraction of ionizing 
photons from a massive halo. We have performed a cosmological RHD simulation,*^ 
where the radiative transfer of UV radiation from individual stellar particles is solved 
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together with the hydrodynamics. In order to evaluate an ionizing photon produc- 
tion rate that depends on stehar age, we have calculated population synthesis with 
PEGASEP^^ in advance of the simulation, and used the obtained rate in the simula- 
tion. The masses of the SPH and dark matter particles are respectively 4.8 x W^Mq 
and 2.9 X 10^ Mq. Thus, we can resolve low mass halos down to M^aio ^ 3 x lO^M© 
with 100 particles. For the reference case, we have also performed a pure hydrody- 
namic (non-RHD) simulation for the same initial conditions and cosmological pa- 
rameters, where the UV radiation transfer is not coupled and therefore no effect by 
UV feedback is included. In a non-RHD simulation, the escape fraction of ionizing 
photons is assessed by solving the UV radiative transfer as a post process using the 
density distributions obtained by a pure hydrodynamic simulation. 

We choose the most massive halo at z = 6.0 in each of the simulations. Hereafter, 
we call the halo in RHD simulation Halo-R, and that in the non-RHD simulation 
Halo-H. The properties of these halos are summarized in Table HITl Here, we should 
mention that the halos are sufficiently massive and hardly lose their baryonic com- 
ponents via photo-evaporation. Actually, the resultant ratio of the baryon mass to 
the halo mass is M^/Mhaio ^ 0.16, which is comparable to Qi^jQ^ = 0.18. 

Fig. [M] shows the numerical results, where number densities, ionized fractions, 
and H2 fractions are shown. For Halo-R, the ionized structures are calculated in 
the RHD simulation, and the escape fractions of HI, Hel, and Hell ionizing photons 
are obtained. On the other hand, in order to assess the escape fractions of ioniz- 
ing photons for Halo-H, we solve the radiative transfer as post-processing until the 
ionization equilibrium is reached at all positions. Hereafter, we call this procedure 
the post-proccesing RT. This is a method frequently used to estimate the escape 
fraction.SSJ'HB \Ye count the number of escaping photons at 2 x 64 x 64 virtual grid 
points placed outside the halo. The upper panels in Fig. [14] show the RHD simula- 
tion, while lower panels are the non-RHD simulations with the post-processing RT 
calculation. 

As seen in Fig. [HI a high density peak forms near the center in Halo-R, where the 
gas is self-shielded from ionizing radiation and H2 molecules form efficiently. This 
central peak allows further star formation. The gas in remaining regions is quite 
smoothly distributed and mostly ionized. This is thought to be the consequence 
of radiation hydrodynamic feedback by phot o-ionizat ion of H, photo-dissociation of 
H2, and photo-heating. On the other hand, several density peaks appear in Halo-H, 
where the gas in each peak is self-shielded from ionizing radiation. Recently, Fernan- 



Table III. Properties of simulated halos 



Name 


Mhalo[M0] 


Mb [Mo] 


M*[M0] 


iV,nt[l/s] 


^e,HI 


Je,Hel 


Je,HeU 


Halo-R (w/ UV) 
Halo-H (w/o UV) 


6.8 X 10^ 
6.7 X 10^ 


1.2 X 10^ 
1.0 X 10^ 


2.9 X 10^ 
1.0 X 10^ 


5.9 X 10^^ 
3.6 X 10^^ 


0.31 

0.15 


0.35 

0.14 


0.13 
0.11 



a) Escape fraction of photons ionizing HI. 

b) Escape fraction of photons ionizing Hel. 

c) Escape fraction of photons ionizing Hell. 
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Fig. 14. Maps of hydrogen number densities, HI number densities and H2 fractions are shown from 
left to right. The shce thickness is 2kpc. The upper panels show the the halo taken from the 
RHD cosmological simulation (Halo-R), while the bottom panels do the halo taken from the 
non-RHD cosmological simulation (Halo-H) with the post-processing RT calculation done by 
the procedure described in the text. 



dez & Shull (2011) have pointed out that the escape fraction sensitively depends on 
the internal structure in a galaxy.^S The clumpy structure in the present simulation 
can reduce the escape fraction of ionizing photons from the galaxy. The resultant 
escape fractions averaged over all directions are presented in Table |TTll The escape 
fractions from the Halo-R are actually higher than those from the Halo-H by a factor 
of 2. It has been often argued that UV feedback would increase the escape fractions 
for a less massive halo because of the mass loss by UV heating.HS Interestingly, the 
present results show that the escape fractions can be enhanced by UV feedback even 
in a massive halo. As for the reionization, what is important is the product of the 
intrinsic ionizing photon production rate TVint multiplied by the escape fraction f^. 
The product is three times higher in the Halo-H, since the star formation rate is 
significantly higher in the Halo-H. 

The dependence of UV feedback on the galaxy mass is unclear in the present 
analysis, and hence a statistical study is needed. But, naively the impacts of UV 
feedback on the escape fractions and star formation are expected to be more dra- 
matic for less massive haloes. Therefore, the present results demonstrate that RHD 
simulations are requisite to reveal the sources of the cosmic reionization. 
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§7. Summary 

In this paper, we have explored the formation of first objects and first galaxies 
with radiative feedbacks. The mass of first objects has been investigated by high- 
resolution hydrodynamic simulations. As a result, it has been revealed that the dark 
matter cusp potential, if it is resolved well, allows a smaller mass object to collapse, 
and then the first object mass could be reduced by roughly two orders of magnitude 
than in previous predictions. This implies that the number of Pop III stars formed 
in the early universe could be significantly larger than hitherto thought. Further- 
more, radiation hydrodynamic simulations with the feedback by photo-ionization 
and photo-dissociation photons have shown that multiple stars can form in a 10^ Mq 
halo. Also, the analytic criteria for radiation hydrodynamic feedback are presented. 

Besides, the fragmentation of an accreting disk around a first protostar has been 
studied with photo-dissociation feedback. As a result, it is found that the frag- 
mentation is not quenched by photo-dissociation feedback, but the "heated bubble" 
resulting from the photo-dissociation can reduce the mass accretion rate onto pro- 
tostars. Also, protostars as small as O.SMq can be ejected through gravitational 
three-body interaction and evolve with keeping their mass. Such small stars might 
be detected as "real first stars" in the Galactic halo. 

Finally, radiation hydrodynamic simulations have been performed to assess the 
impact of UV feedback on first galaxies. By comparing radiation hydrodynamic 
simulations to pure hydrodynamic simulations, we have found that UV feedback 
enhances the escape fraction by a factor of 2 even in a massive halo. But, the 
product of the intrinsic ionizing photon production rate multiplied by the escape 
fraction is three times higher in a case without UV feedback than in a case with UV 
feedback. These results imply that UV feedback deserves careful consideration to 
reveal the cosmic reionization history. 
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